

%% load estimation results
load('parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
clearvars -except thetaStar K mpow N reb X Y YY Z analyst specif


%% estimated mean utilities
UStar = payoffs(thetaStar,X,Z,YY,reb,mpow);


%% counterfactuals
load('counterfactualUK_replication.mat','srng_cfactUK')
rng(srng_cfactUS)
NS = 2e+3;
YY4 = YY(YY(:,6)==0,1:5); % 4-mpow cw-intervention strategy profiles
timeUK=0;

% UK never intervenes
VnoUK_R = UStar(YY(:,3)==0,1,1:N) + randn(size(YY4,1),NS,N); % entry payoffs of rebel group
VnoUK_US = UStar(YY(:,3)==0,2,1:N) + randn(size(YY4,1),NS,N);
VnoUK_UK = UStar(YY(:,3)==0,3,1:N) + randn(size(YY4,1),NS,N);
VnoUK_France = UStar(YY(:,3)==0,4,1:N) + randn(size(YY4,1),NS,N);
VnoUK_Russia = UStar(YY(:,3)==0,5,1:N) + randn(size(YY4,1),NS,N);
VnoUK_China = UStar(YY(:,3)==0,6,1:N) + randn(size(YY4,1),NS,N);
filen=1:NS;
noUKintervEq=cell(NS,N);
for n=1:N
    tic
    parfor b=1:NS
        fn=filen(b);
        fname=['gameUK',num2str(fn),'.nfg'];
        fileID=fopen(fname,'w');
        fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
        fprintf(fileID,'{ "US" "France" "Russia" "China" } { 2 2 2 2 }\n',1); %#ok<CTPCT>
        formatSpec='%12g ';
        fprintf(fileID,formatSpec,reshape([VnoUK_US(2:end,b,n),VnoUK_France(2:end,b,n),VnoUK_Russia(2:end,b,n),VnoUK_China(2:end,b,n)]',1,4*2^4)); %#ok<*PFBNS>
        fclose(fileID);        
        [~,NE]=system(['gambit-enumpoly ',fname]);
        findNE=strfind(NE,'NE');
        EQ=zeros(2^4,length(findNE));
        EV_R=zeros(1,length(findNE));
        for m=1:length(findNE)
        	sigma_US=str2double(NE(findNE(m)+12:findNE(m)+19));
            sigma_France=str2double(NE(findNE(m)+30:findNE(m)+37));
            sigma_Russia=str2double(NE(findNE(m)+48:findNE(m)+55));
            sigma_China=str2double(NE(findNE(m)+66:findNE(m)+73));
            EQ(:,m)=(sigma_US.^YY4(2:end,2)).*((1-sigma_US).^(1-YY4(2:end,2))).*...
                (sigma_France.^YY4(2:end,3)).*((1-sigma_France).^(1-YY4(2:end,3))).*...
                (sigma_Russia.^YY4(2:end,4)).*((1-sigma_Russia).^(1-YY4(2:end,4))).*...
                (sigma_China.^YY4(2:end,5)).*((1-sigma_China).^(1-YY4(2:end,5)));
            EV_R(m)=VnoUK_R(2:end,b,n)'*EQ(:,m);
        end
        noUKintervEq{b,n}=[EV_R;EQ];
    end
	toc
    timeUK=timeUK+toc;
end
[noUKDMES,noUKintervEq,noUKeqPayoffs] = createDMES(noUKintervEq,VnoUK_R,VnoUK_US,VnoUK_UK,VnoUK_France,VnoUK_Russia,VnoUK_China,[YY4(:,1:2),zeros(size(YY4,1),1),YY4(:,3:5)],specif);

% UK always intervenes
VallUK_R = UStar([1;find(YY(:,3)==1)],1,1:N) + randn(size(YY4,1),NS,N); % entry payoffs of rebel group
VallUK_US = UStar([1;find(YY(:,3)==1)],2,1:N) + randn(size(YY4,1),NS,N);
VallUK_UK = UStar([1;find(YY(:,3)==1)],3,1:N) + randn(size(YY4,1),NS,N);
VallUK_France = UStar([1;find(YY(:,3)==1)],4,1:N) + randn(size(YY4,1),NS,N);
VallUK_Russia = UStar([1;find(YY(:,3)==1)],5,1:N) + randn(size(YY4,1),NS,N);
VallUK_China = UStar([1;find(YY(:,3)==1)],6,1:N) + randn(size(YY4,1),NS,N);
filen=1:NS;
allUKintervEq=cell(NS,N);
for n=1:N
    tic
    parfor b=1:NS
        fn=filen(b);
        fname=['gameUK',num2str(fn),'.nfg'];
        fileID=fopen(fname,'w');
        fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
        fprintf(fileID,'{ "US" "France" "Russia" "China" } { 2 2 2 2 }\n',1); %#ok<CTPCT>
        formatSpec='%12g ';
        fprintf(fileID,formatSpec,reshape([VallUK_US(2:end,b,n),VallUK_France(2:end,b,n),VallUK_Russia(2:end,b,n),VallUK_China(2:end,b,n)]',1,4*2^4)); %#ok<*PFBNS>
        fclose(fileID);        
        [~,NE]=system(['gambit-enumpoly ',fname]);
        findNE=strfind(NE,'NE');
        EQ=zeros(2^4,length(findNE));
        EV_R=zeros(1,length(findNE));
        for m=1:length(findNE)
        	sigma_US=str2double(NE(findNE(m)+12:findNE(m)+19));
            sigma_France=str2double(NE(findNE(m)+30:findNE(m)+37));
            sigma_Russia=str2double(NE(findNE(m)+48:findNE(m)+55));
            sigma_China=str2double(NE(findNE(m)+66:findNE(m)+73));
            EQ(:,m)=(sigma_US.^YY4(2:end,2)).*((1-sigma_US).^(1-YY4(2:end,2))).*...
                (sigma_France.^YY4(2:end,3)).*((1-sigma_France).^(1-YY4(2:end,3))).*...
                (sigma_Russia.^YY4(2:end,4)).*((1-sigma_Russia).^(1-YY4(2:end,4))).*...
                (sigma_China.^YY4(2:end,5)).*((1-sigma_China).^(1-YY4(2:end,5)));
            EV_R(m)=VallUK_R(2:end,b,n)'*EQ(:,m);
        end
        allUKintervEq{b,n}=[EV_R;EQ];
    end
	toc
    timeUK=timeUK+toc;
end
[allUKDMES,allUKintervEq,allUKeqPayoffs] = createDMES(allUKintervEq,VallUK_R,VallUK_US,VallUK_UK,VallUK_France,VallUK_Russia,VallUK_China,[YY4(:,1:2),ones(size(YY4,1),1),YY4(:,3:5)],specif);


%%
clearvars -except N NS srng_cfactUK YY4 timeUK VnoUK_R VnoUK_US VnoUK_UK VnoUK_France VnoUK_Russia VnoUK_China ...
     VallUK_R VallUK_US VallUK_UK VallUK_France VallUK_Russia VallUK_China noUKintervEq noUKDMES noUKeqPayoffs allUKintervEq allUKDMES allUKeqPayoffs
save('counterfactualUK.mat')



